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We report on our experience with the two-step multi-bosonic algorithm in a large scale Monte Carlo simulation 
of the SU(2) Yang-Mills theory with dynamical gluinos. First results are described on the low lying spectrum of 
bound states, the string tension and the gluino condensate. 


1. INTRODUCTION 

An important step towards the understanding 
of non-perturbative properties of supersymmet¬ 
ric (SUSY) gauge theories is the numerical simu¬ 
lation of the supersymmetric extension of Yang- 
Mills (SYM) theory. The theoretical motivation 
for this has been summarized at the previous [Q 
and present [f| lattice conferences. (For refer¬ 
ences see these reviews and ||].) First steps to¬ 
wards a numerical Monte Carlo simulation with 
dynamical gluinos have already been presented a 
year ago [Q. Since then our collaboration made 
important progress in a first large scale simulation 
of SU(2) SYM on the CRAY T3E-512 at HLRZ 
Jiilich. This is a short status report which will be 
followed soon by a more detailed publication. 

1.1. Lattice action 

The SU(2) Yang-Mills theory with gluinos (= 
Majorana fermions in the triplet representation) 
becomes supersymmetric in the massless limit. 
Massive gluinos break SUSY softly. For a lattice 
regularization one can take the Wilson action for 
gauge field and gluinos, as proposed some time 
ago by Curci and Veneziano ||. This contains 
two bare parameters: (3 for the gauge coupling 
and K for the hopping parameter (bare gluino 


mass). The Majorana nature of the gluino is 
taken into account by considering Nf = tj ad¬ 
joint Dirac flavours. The effective action for the 
gauge field assumes the form 

Scv = - \ TrU pi) ~ \ l°g det Q[U] ,(1) 

P i ' ' 

where the fermion matrix is 


Q 


yv,xu — 3y T 8 v 


K X + 'Y [i)Vvu ,Xfl (2) 

m=± 


with the gauge link in the adjoint representation 
Vvu,xp — 2 Yr (U^.^T V U x pTu) • 

An interesting new development is that an ex¬ 
actly zero mass gaugino can be described on the 
lattice by the Neuberger- action Q. (Concerning 
chiral symmetry see 0.) The technical difficulty 
is to determine the necessary inverse square-root 
(Qt(5) _ 2. For this we can use the quadratically 
optimized polynomials discussed in sec. |j. The 
results of a numerical study of the Neuberger- 
action will be published elsewhere. 


1.2. Pfaffians 

The Curci-Veneziano action assumes det(Q) 3 
for the Majorana fermion. This may lead to a sign 
problem because the path integral for Majorana 
fermions gives the Pfaffian 

Pf(M) = [ 
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jyi'iN^ ol i0i--- ol nPnM ai f3i ■ ■ ■ M aN p N , (3) 

where M = CQ satisfies 

Pf(M ) 2 = det(Af) = det(Q) = det(Q) . (4) 

Here Q = 75 Q is the hermitean fermion matrix. 
Since Q has doubly degenerate real eigenvalues, 
the determinant is non-negative: det(Q) > 0 . 
The sign of the Pfaffian is, however, unknown. 

A numerical procedure for the computation of 
Pfaffians can be based on the decomposition || 

M = P T JP 1 Pf(M) = det(P) , (5) 

where J is a block-diagonal matrix containing on 
the diagonal 2 ® 2 blocks equal to e = * 02 . This 
form of M can be achieved by a procedure analo¬ 
gous to the Gram-Schmidt orthogonalization and 
then P turns out to be a triangular matrix. (See 
e.g. the treatment of simplectic groups in |j|.) 

The numerical procedure requires the storage 
of a full N (g) N matrix, similarly to the deter¬ 
minant calculation by LU-decomposition. There¬ 
fore, one can only deal with relatively small lat¬ 
tices. In a test with dynamical gluino updating 
at (/3 = 2.3, K = 0.1925) on 4 3 • 8 lattice the 
Pfaffian of every randomly chosen configuration 
turned out to be positive (see fig. jl"). 
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Figure 1. The values of the pfaffian on 4 3 • 8 
lattice at (/? = 2.3, I< = 0.1925) versus the largest 
eigenvalue \ max of Q^Q. 


In general, however, the Pfaffian is not al¬ 
ways positive. In fact, if a pair of degenerate 
eigenvalues of the hermitean fermion matrix Q 
changes sign, it is plausible that the sign of Pf (M) 
changes, too. An example is shown by fig. gj, 
where one of the configurations in fig. [l] is con¬ 
sidered as a function of the “valence” hopping 
parameter K v in M. As the numerical determi¬ 
nation of the nearly zero eigenvalues of Q shows, 
at the same K v where Pf(M) changes sign there 
is also a sign change of an eigenvalue pair. 


lattice: 4®4®4®8 



K v 

Figure 2. The absolute value of the pfaffian on 
a 4 3 • 8 configuration as a function of the hopping 
parameter in M (and Q ). Open points stand for 
Pf(M) > 0, full ones for Pf(M) < 0. 

The spectral flow of Q is relevant for the 
overlap-inspired fermionic definition of the topo¬ 
logical charge (see |To| and references therein). 
Since the index of a massless Dirac operator in the 
adjoint representation of SU(2) in a gauge field 
background of topological charge Q top is 4 Q top 
p| , this suggests that at the hopping parameter 
value K' used in the Neuberger action the sign of 
Pf(M) is given by 

Pf(M(K'))/\Pf(M(K'))\ = e 2niQt °r. ( 6 ) 

One has to have in mind that the fermionic defini¬ 
tion can also give half-integer topological charges. 
(For SU(Ac) gauge theory the topological charge 
can be an integer multiple of l/N c .) Since K' is 
larger than the dynamical hopping parameter K 
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and the level crossings typically occur between K 
and K', it is plausible that the path integral at 
K is dominated by configurations with positive 
Pfaffian. Therefore the Curci-Veneziano action 
gives the same continuum limit as the definition 
with the path integral over Majorana fermions. 

The question of fractional topological charges is 
important for the low energy effective action [[L3| 
which assumes a remnant chiral symmetry Z 2 n c - 
In fact, if fractional topological charges would ex¬ 
ist, the ^-parameter would be periodic only with 
27 tN c and the remnant chiral symmetry would be 
Z 2 , instead of Z 2 n c ■ In this case the spontaneous 
symmetry breaking Z 2 n c —> Z 2 would be absent 
and there would be no first order phase transition 
at zero gluino mass. In the quenched continuum 
limit the half-integer topological charges seem to 
persist p~2| but the continuum limit with dynam¬ 
ical gluinos may be different. 


The difficulty for small fermion masses is that 
the condition number \/e becomes very large 
(10 4 — 10 6 ) and very high orders n = 0(1O 3 ) 
are required for a good approximation. This re¬ 
quires large storage and the autocorrelation be¬ 
comes bad since it is proportional to n. 

2.1. Improved version: two-step updating 

Substantially smaller storage and shorter au¬ 
tocorrelations for small fermion masses can be 
achieved by starting from a two-step approxima¬ 
tion Jl5[ : 

{det(Q t Q)} Ar//2 


" detP^ ) (QtQ)detP^ ) (QtQ) 

where now 

lim P«(x)P^0r) = x~ N t' 2 , x€[e, A] .(12) 

n 2—>00 


2. MULTI-BOSONIC ALGORITHM 

The multi-bosonic algorithm for Monte Carlo 
simulations of fermions has been proposed by 
Liischer 0- In the original version for Nf 
flavours one considers the approximation of the 
fermion determinant 

|det(Q)r / = (det^Q )}^ 72 
1 

detP„(Q t Q) 

where the polynomial P n satisfies 
lim P n (x) = x~ N f /2 

n —>oo 

in an interval [e, A] covering the spectrum of Q^Q. 
For the multi-bosonic representation of the deter¬ 
minant one uses 

n 

P n (QtQ) = P„ ( Q2 } = ro J](Q _ p *)(Q _ p .) (9) 


JJdet [(Q - Pj)(Q - pj)] 1 oc 
3 =1 


(7) 

( 8 ) 


[d<f>] exp{— 


3 =1 X V 


The multi-bosonic representation of the determi¬ 
nant is used for det P},)^. The correction factor 
det Pn? is realized in a noisy correction step |l6| 
with the accept-reject function 

where rj is generated from the simple Gaussian 
noise rj with a suitable polynomial approximation 
Pnf as 

77 = P${Q[U?Y*rf =► p£HQ[U\ 2 W ■ (14) 

An important gain in performance is obtained 
by the use of quadratically (“least-square”) op¬ 
timized polynomials [jl7|]. The Chebyshev- (for 
Nf = 2) or Legendre- (for Nf = 1) polynomials 
are bad for large A/e. The quadratically opti¬ 
mized polynomials a re m uch better. Note that, 
as discussed in sec. u for gluinos one has to 
consider Nf = 

A complete cycle of sweeps contains heatbath 
and overrelaxations sweeps for the boson fields 
and a Metropolis sweep for the gauge field fol¬ 
lowed by the accept-reject step. We choose the 
order of the first polynomial such that the ac¬ 
ceptance rate in the accept-reject step is typi¬ 
cally 80-90%. The longest autocorrelations ap¬ 
pear for “gluonic” quantities depending on gauge 
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field variables. For instance, in our runs the auto¬ 
correlation of the plaquette expectation value is 
typically given by ~ 300 — 400 cycles. The 
“fermionic” autocorrelations are about a factor 
~ 10 smaller. 

2.2. Measurement correction 

To avoid systematic errors, both in the original 
version and in the two-step version the polyno¬ 
mial approximations have to be taken, in prin¬ 
ciple, to infinite order. In praxis it is enough 
that the errors introduced by the polynomial 
approximations are smaller than the statistical 
errors. However, due to the limited approxi¬ 
mation quality, for small fermion masses some¬ 
times exceptional configurations with very small 
(1CT 8 — 10” 10 ) eigenvalues of Q^Q are generated. 
In order to avoid these one could increase the 
polynomial orders but it is more economical to 
perform a measurement correction step by a gen¬ 
eralization of the method of ref. [[L 8 |. In this way 
the expectation value of a quantity A is given by 

= (Aexp{r?t[l - , 15 , 

(exp { 77 ! [1 - PnJ (Q^Q)}v})u,ri 

where rj is a Gaussian noise and the polynomial 

is such that P$ P^ P$ is an optimized ap¬ 
proximation of the function x~ Nf / 2 in the inter¬ 
val [e',A]. This correction step can be arranged 
in such a way that one obtains an “exact” algo¬ 
rithm. For instance, e' can be set to zero and then 
the correction step is not sensitive to the fluctu¬ 
ations of the smallest eigenvalue. (Note that in 
our case the least-square optimization can also 
be achieved for zero lower boundary of the inter¬ 
val, in contrast e.g. to the minimization of the 
maximal relative deviation with Chebyshev poly¬ 
nomials.) For the evaluation of P$ one can use 
714 -independent recursive relations 0 which can 
be stopped by observing the convergence of the 
result. Another possibility is to take e' > 0 and to 
perform the correction exactly for the few excep¬ 
tional eigenvalues below e' by multiplying with 
\ N f/ 2 P$ (\)pf$ (A). 

Experience shows that for most physical quan¬ 
tities the exceptional configurations give no ex¬ 
ceptional contributions. There are, however, also 


some quantities which are sensitive to the small 
eigenvalues, hence the measurement correction is 
essential (see fig. ||). As the figure shows, even 


Lattice: 6*3 x 12, K=0.195, beta=2.3 

_ corrected data 




raw data 



Figure 3. The measurement correction for the a- 
pion propagator at zero distance. The exceptional 
configurations with small eigenvalues contribute 
strongly to the raw data. After correction these 
contributions are still important but of normal 
size. 


the small lower limit of the polynomial approxi¬ 
mation interval e = 10 5 , in connection with the 
polynomial orders ni = 24, 712 = 200, is insuf¬ 
ficient for suppressing the exceptional configura¬ 
tions. In this case a measurement correction step 
with e' = 0 and 714 = 400 is enough for good 
precision. 


3. APPROACHING 
SUPERSYMMETRY 

The Monte Carlo simulations have been per¬ 
formed up to now mainly at P = 2.3. After 
last years tests at lower hopping parameter val¬ 
ues [0 we are presently concentrating on the in¬ 
terval 0.185 < K < 0.1975. Lattice sizes range 
from 4 3 • 8 to 12 3 • 24. The characteristic fea¬ 
ture of these simulations are the rather large con¬ 
dition numbers growing up to A/e = 10 5 — 10 6 . 
Our present best estimate for the critical hopping 
parameter for zero gluino mass is Kq — 0.195 
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(see sec. |ib||). This means that the bare gluino 
mass in lattice units at e.g. K = 0.1925 is about 
amo = ^(Kq 1 — A" -1 ) ~ 0.03. This gives a ra¬ 
tio to the smallest glueball mass mo/M 0+ ~ 0.08. 
The large condition numbers are further increased 
if the “valence” hopping parameter K v in the 
gluino propagators is increased. A typical ex¬ 
ample for the average of the lowest eigenvalue is 
shown in fig. [|. (The approximation interval for 


lattice: 8®8<8>8<8>16 



Figure 4. The average lowest eigenvalue of Q^Q 
as a function of the inverse (valence) hopping pa¬ 
rameter in a dynamical update at (/3 = 2.3, K = 
0.1925) on 8 3 • 16 lattice. The error bars here show 
standard deviations (not the errors of averages). 

polynomials is in this case [e, A] = [0.0001,3.7].) 
As one can see, near 1 /K v = 1/0.2025 ~ 4.94 
the average smallest eigenvalue becomes rather 
small and fluctuates practically down to zero. For 
smaller values of 1/K V the picture remains un¬ 
changed. This agrees with the findings of ref. [|Io| 
about the absence of the spectrum gap for a range 
of bare fermion masses. (Note that the start of 
the gapless region at K v = 0.2025 lies higher than 
K 0 because in fig. |l| the dynamical “sea” gluino 
mass is not changing.) 

3.1. Bound state masses 

A basic assumption about the dynamics of 
SYM theory is that there is confinement. The 
colourless bound states should be organized in 
degenerate supermultiplets which are split up at 
non-zero gluino mass. The pattern of states 


should be similar to QCD with a single flavour of 
quarks. For instance, the low-lying meson states 
have quantum numbers as rf, fo etc. Since our 
constituent gluinos are in the adjoint representa¬ 
tion, we call these states a-rf, a-fo etc. There 
are, of course, also the glueball states ( gg) as in 
pure gauge theory and for completing the super¬ 
multiplets we also need spin-^ states: the gluino- 
glueballs (gg). 

We determine the masses of low-lying states 
by using the blocking-smearing technology for 
achieving better projections to the sources and 
sinks. Our preliminary results are summarized 
in fig. ||. The presence of two distinct light 0 + 



Figure 5. The preliminary results for bound 
state masses. 

states in the spectrum seems to prefer the gen¬ 
eralization of the Veneziano-Yankielowicz action 
@ proposed in ref. |jl9f . 

3.2. String tension 

The confinement can be characterized by the 
string tension obtained from the linear part of the 
static potential at large distances. We determined 
the static potential from APE-smeared Wilson 
loops following m . At (/? = 2.3, I\ = 0.19) on 
8 3 • 16 lattice the optimized smearing radius is 
3.3 and we obtained stable results against varia¬ 
tions of the T fit-interval. As shown by fig. |[ the 
result for the square root of the string tension in 
lattice units is a^fo = 0.21(1). At K = 0.1925 the 
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Figure 6. The three-parameter fit to the static 
potential: V(R) = Vo + aR — a/R at (/? = 

2.3, K = 0.19) on 8 3 • 16 lattice. 

corresponding value is surprisingly small: ay/a = 
0.10(3) but the 8 3 • 16 lattice might be too small 
in this case, causing finite volume effects. 

3.3. Gluino condensate 

If there are no fractional topological charges, 
the hopping parameter Kq corresponding to zero 
gluino mass is signaled by a first order phase tran¬ 
sition which is due to the spontaneous discrete 
chiral symmetry breaking Z 4 —> Z 2 . The gluino 
condensate should have a jump of order 0 (1) in 
lattice units ||. It is possible that this transition 
develops only in the continuum limit and at finite 
lattice spacings it shows up just as a fast cross¬ 
over. Our preliminary results are summarized in 
fig. [t| The jump at 1/K = 1/0.195 ~ 5.13 sug¬ 
gests a vanishing gluino mass at Kq ~ 0.195. 
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sqrt(ff) = 0.2201 (121), [1,4]/ 
sqrt(ff) = 0.20 6 (118), [ ,5] ( 
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